The Two-Dimensional Disordered Boson Hubbard Model: Evidence for a Direct 
Mott Insulator-to-Superfluid Transition and Localization in the Bose Glass Phase 
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We investigate the Bose glass phase and the insulator-to-superfluid transition in the two-dimensional 
disordered boson Hubbard model in the Villain representation via Monte Carlo simulations. In the 
Bose glass phase the probability distribution of the local susceptibility is found to have a 1/x 2 tail 
and the imaginary time Green's function decays algebraically C(r) ~ r _1 , giving rise to a divergent 
global susceptibility. By considering the participation ratio it is shown that the excitations in the 
Bose glass phase are fully localized and a scaling law is established. For commensurate boson 
densities we find a direct Mott insulator to superfluid transition without an intervening Bose glass 
phase for weak disorder. For this transition we obtain the critical exponents z = l,v = 0.7 ± 0.1 
and r\ = 0.1 ± 0.1, which agree with those for the classical three-dimensional XY model without 
disorder. This indicates that disorder is irrelevant at the tip of the Mott-lobes and that here the 
inequality v > 2/d is violated. 



PACS numbers: 67.40-w, 74.70.Mq, 74.20.Mn 
I. INTRODUCTION 

The localization transition in disordered systems of 
interacting bosons has attracted a lot of interest in re- 
cent years. The observation is that there is a direct 
superconductor to insulator transition at zero temper- 
ature, which can be tuned by varying a suitable con- 
trol parameter like the disorder strength or an applied 
magnetic field. The paradigmatic two-dimensional real- 
izations of such systems are amorphous superconducting 
filmsn, where the Cooper pairs are considered as bosons, 
and 4 He adsorbed in porous median On the theoretical 
side, these systems are described most simply by a boson 
Hubbard model, and a large amount of effort has been 
put forth to study the various aspects of this transition. 
This includes mean-field xalculationgj, real space renoija 
malization group studiesu, strong coupling expansionsu 
and quantum Monte Carlo investigations in one-El and 
two dimensionsQij. 

Of particular interest are the properties of the so called 
Bose glass phase which is supposed to emerge for strong 
disorder, and which is characterized by being insulating 
but compressible and gaplesstl The latter property leads 
to a diverging superfluid susceptibility throughout the 
Bose glass phase, which is reminiscent of the quantum 
Griffiths phase in random transverse Ising systemsliHTL-X 
It is one of the main goals of this work to investigate the 
Bose glass phase and to characterize the Griffiths singu- 
larities by means of Monte Carlo simulations of the 2D 
Bose Hubbard model. The other major point of interest 
is the insulator to superfluid transition for weak disor- 
der and commensurate boson densities, since there is an 
ongoing debate whether this transition occurs fro 
Mott insulating (MI)- or the Bose glass (BG) phas< 
In this paper we present numerical evidence for a direct 



Mott insulator to superfluid transition that falls in the 
universality class of the pure classical XY-model in three 
dimensipas. A short account of our results has been given 
recently^ . 

The paper is organized as follows: In Sec. [n] the 2D 
Bose Hubbard model is introduced and it is briefly stated, 
how one arrives via a set of transformations at a 3D classi- 
cal model, which is used subsequently in the simulations. 
In Sec. lit we introduce the quantities of interest and 
summarize the scaling theory of Ref.0 which is needed 
for a finite size scaling analysis of the data. In Sec. |y| 
the Monte Carlo algorithm is presented. In Sec. [v] we 
examine the Bose glass phase for strong disorder and in 
Sec. VI we consider the case of weak disorder and address 



the question of a direct MI-SF transition. 



II. THE MODEL 

We consider the boson Hubbard model with disorder 
on a square lattice in d = 2 dimensions. The Hamiltonian 
is given by 



Hbh = Hq + Hi 



H = - t J2(&A + <Mj) 

(id) 



Hi = — ^n? - + v i)^ > 



(1) 



(2) 



(3) 



where $i ($|) is the boson annihilation (creation) opera- 
tor at site i and hi = is the usual number operator. 
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The hopping matrix is t for nearest neighbors and zero 
otherwise. The bosons interact via an on-site repulsion 
U, fi is the chemical potential and represents the ran- 
dom on-site potential varying uniformly in space between 
-A and A. 

As a simplification we take the "phase-only" 
approximation!!], where amplitude fluctuations in (Q) are 
integrated out. One assumes that the complex field <f> has 

the simple form $i = |$o| e i where fa is the phase op- 
erator conjugate to the number operator hi with commu- 



tation relation 



phase-Hamiltonian 



] = iSij. This yields the quantum- 



Hqph = H + H i 



cos 



{id) 



Hi = IT ™i ~ + Vl 



(4) 



(5) 



(6) 



where hi now measures the deviation from a mean den- 
sity. 



A. Duality transformation 

Through a sequence of transformations the 2D quan- 
tum mechanical Hamiltonian can be mapped onto a 
classical model of divergence-free integer link variables 
in (2 + f ) dimensions which is suitable for Monte Carlo 
simulations!! _ . 

The transformation uses the standard Suzuki- TrottertZI 
formula where the inverse temperature (3 is divided into 
L r time slices of width At = f3/L T yielding a path inte- 
gral representationj-af Hamiltonian (|J). Taking the Vil- 
lain approximation!!!!, which replaces the exponentiated 
cosine term in the partition function by a sum of peri- 
odic gaussians, and after a further Poisson summation, 
one can show that the ground state energy density of the 
quantum model (^) is equal to the free energy density of 
a purely classical model 



/ 



T 
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-PHqph 
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L 2 L t At 



where the classical Hamiltonian i? c iass is given by 

VJ=0 ,^ > 

ffclass = I 2 J ^ T ' ~ ^ + V *) J *,-r' 

(i,r') V 



(7) 



(8) 



and the parameter K acts as an effective temperature of 
the model and is used to tune through the phase transi- 
tion. K corresponds to t/U in the quantum phase model 



(Q) and thus by changing K one moves horizontally, i.e. 
parallel to the t/U axis in the phase diagram of Hamilto- 
nian (^). The J's = (J x , J y , J T ) are three-component 
vectors consisting of integer variables on the links of 
the original lattice running from — oo to oo. The r- 
component J T corresponds to the particle number. Note 
that the summation in (||) is just over divergence-free 
configurations of the link variables. Strictly speaking, 
one has to take the limit L T — > oo for the equality in 
eq. (^) to be valid. In the simulations however, we have 
to keep L T finite, which means that we are working at 
a temperature T ~ 1 / L T for the original quantum phase 
model (§. 

The action (||) describes an ensemble of string like 
objects in a random potential, and considering the 
divergence-free condition, these objects can be inter- 
preted as flux-lines in a high-Tc superconductor. This 
observation is not surprising, since one has the general 
correspondence between 2D bosons at zero temperature 
and flux-lines in 3D high- Tc superconductorscjE!!. 



B. Phase diagram 

The phase diagram of the classical Hamiltonian (||) in 
the thermodynamic limit (L — > oo, L T — > oo), which 
corresponds to the T = phase diagram of the quan- 
tum model (Q) can be obtained qualitatively by simple 
perturbative argumentsu. 

For the case without disorder (Vi — 0) one obtains 
two phases: a Mott insulating and a superconducting 
phase. In the /i vs. K plane the Mott insulating phase 
has a lobe like shape which can be understood by first 
considering the atomic limit K — and then allow for 
small hopping in the quantum phase model (0), i.e. finite 
effective temperature for the classical model (g). At zero 
effective temperature, every site is occupied by an integer 
number n = J T of bosons which minimizes the on site 
energy e(n) = —/in + \v? . Thus in the interval n < 
/i+5 < n+1 the occupation number is fixed to n at every 
site causing the compressibility k — dn/d^i to vanish. 
Furthermore there is a gap in the single particle spectrum 
making the Mott phase indeed an insulator. 

For increasing hopping amplitude t/U the bosons in 
the original quantum phase model can gain kinetic 
energy by hopping on the lattice and thus the gap de- 
creases with increasing t. When the gap vanishes, the 
bosons are free to hop on the lattice thus producing su- 
perfluidity. Equivalently in the classical model (g) an 
increase in effective temperature K causes a distortion of 
the initially straight world lines allowing them to move 
around freely above some critical value for K giving the 
Mott phase a lobe like shape. As Hamiltonian (||) is pe- 
riodic in /j, with period f , the lobes are repeated on the 
fi axis. 

In the presence of disorder a-pew, insulating Bose glass 
phase is supposed to emergeEl in addition to the Mott 
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insulating and the superfluid phase. As the disorder 
strength A increases, the Mott lobes shrink since there 
will always be sites where the energy necessary to add a 
particle or hole is reduced due to the disorder and these 
sites therefore favour the addition of particles or holes. 
Thus the compressibility ceases to vanish and the trace 
of the Mott insulator disappears. Fisher et al. now ar- 
gue that the extra bosons will not immediately produce 
superfluidity but will still be localized due to the dis- 
order. Thus the new phase is characterized by a non- 
vanishing compressibility, no energy gap and zero con- 
ductivity. Only for larger hopping amplitudes the bosons 
will become delocalized resulting in a superfluid phase. 

The phase diagram for weak disorder suggested by our 
sinnJations is depicted in Fig. |l|. Note that contrary to 
Ref.cl we find a direct MFSF transition for integer values 
of the chemical potential (see Sec. VI) and therefore 
there is no intervening BG phase at the tip of the lobe. 




superfluid 



superfluid 



K 



FIG. 1. The L T = oo phase diagram of the classical model 
(^) for weak disorder suggested by our results presented in this 
paper. The Mott lobes are centered around integer values of 
the chemical potential. For integer values of fj, we find a direct 
MFSF transition at a multicritical point, therefore there is no 
intervening BG phase at the tip of lobe. Since the Hamilto- 
nian (^) is periodic in the chemical potential with period one, 
the Mott lobes are repeated in the y-direction. The dots de- 
note the phase transitions we examine in the simulations and 
are labeled for later reference. The arrows indicate whether 
/i or K were used to tune through the different transitions. 



disappeared and only the BG- and the superfluid phase 
remain. 
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FIG. 2. The zero temperature phase diagram of the classi- 
cal model Q for disorder strengths A > 1/2. Fhe Mott lobes 
have disappeared and only the BG- and the SF-phase remain. 
The dot denotes the BG-SF transition studied in Sec. [v| and 
is identical to point #3 in Fig. [[]. 

Furthermore Fisher et al. argue that due to the non- 
vanishing density of states at zero energy the imaginary 
time Green's function C(r) decays algebraically in the 
BG phase. This leads to a divergent superfluid suscepti- 
bility x — J C( T ) dr in the BG phase, which thus resem- 
bles the quaii^m—Griffiths phase found in other disor- 
dered systemsErli-i We investigate this scenario in Sec. 



III. SCALING THEORY AND FINITE-SIZE 
SCALING 

For completeness we give a brief summary of the T = 
scaling theory developed in Ref.u. Furthermore we intro- 
duce the quantities of interest and show how they can be 
expressed in terms of the link variables of the classical 
model (||). Since in the simulations we are working with 
finite systems, it is also shown how to obtain results via 
finite-size scaling. 



A. Stiffness 



For disorder strengths A > {7/2 one obtains the phase 
diagram shown in Fig. 0. The Mott lobes have totally 



From the anisotropy of the classical model (|8|) it can 
be seen, that near a continuous phase transition there 
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are two diverging correlation " lengths" £ and £ r . This is 
due to the fact that we are considering a T = quantum 
phase transition where static and dynamic properties of 
a system are linked. £ is the correlation length in the spa- 
tial directions, whereas £ T denotes a diverging time scale 
in the imaginary time direction. As usual we introduce 
the critical exponent v by £ ~ 8~ u . Generally £ r will di- 
verge with a different exponent so we define the dynam- 
ical critical exponent z by £ T ~ £ z . The control parame- 
ter (5 measures the normalized distance from the critical 
point, which is 6 — {K — Kc)/Kc or S = (/x — He) I He 
in the simulations. 

An important quantity in the further analysis is the 
zero frequency stiffness p(0), which is proportional to the 
order parameter, i.e. the condensate. Thus by calculat- 
ing the stiffness we are able to locate the critical point 
of a phase transition from an insulating to a supercon- 
ducting phase. The stiffness measures the change of the 
free energy under imposing a twist 6 in the phase of the 
order parameter and is defined by 



Sfs 



1 



P {vef 



(9) 



where f s denotes the singular part of the free energy. Via 
hyperscaling the scaling form 

p ru £~( d +z-V (10) 

can be derived^. 

To calculate the stiffness in the simulations we have 
to express p in terms of the link variables of the classical 
model (|J). By considering the effect of an external vector 
potential one obtainsu 



P(0) = ^-[(nl)U , 
where n x defines the winding number 

H * = | E J (i, r) > 
(i,r) 



(11) 



(12) 



and [. . .] av denotes a disorder average. Since at the 
critical point the diverging correlation lengths are much 
larger than the system size we need to obtain a finite-size 
scaling form of the stiffness to analyse our data. Because 
there are two diverging correlation lengths, all quantities 
will depend on the ratios L/£ and L T /£ T . By considering 
this and eq. (|o|) we get 



p(0) = r {d+z+2) P(L//;,L T /t T ) 



which can be expressed as 



P(0) = I ^P(L 1/v 6,L r /L*) 



(13) 



(14) 



where P and p are scaling functions. Thus we see that by 
keeping the aspect ratio a = L T /L Z constant, the scaling 



function p just depends on one argument. This enables 
us to locate the critical point, since at criticality 5 = 
and L d+z - 2 p(0) then is independent of L, which means 
that a plot of L d+z ~ 2 p(0) vs. K for different system sizes 
with constant aspect ratio yields an intersection point at 
the critical point. 



B. Compressibility 

In a similar way as for the stiffness, a scaling expression 
for the compressibilty 



dn 
Oft 



(15) 



can be derived. By imposing a twist in the phase of the 
order parameter in the imaginary time direction instead 
of the space directions, one obtains the scaling form 

K~r (d - z) , (i6) 

and the corresponding finite-size scaling form is 



k(0) = -J— k ( lM v 5, — 
K ' L d ~ z V L z 



(17) 



Since by definition n measures the fluctuations in the 
particle number, which is given by the J['s in the clas- 
sical model,-^), k can be expressed in terms of the link 
variables asQ 



k(0) 



1 



L d L T 

where N b = £„ r / J7 T >- 



(18) 



C. Correlation functions 

In addition to v and z a third critical exponent r\ 
can be introduced by considering the correlation func- 
tion C(r, r', r, r') for creating a boson at point (r, r) and 
destructing it at (r',r'). In the quantum phase model 
(0) this correlation function is given by 



C{v,v',t,t') = [(e^'M-kWl)] 



(19) 



and due to the disorder average C(t,t',t,t') will 
be translationally invariant, i.e. C(r, r',r, t 1 ) = 
C(r — r', r — r'). We will only be interested in the equal 
time correlation function C x (r) = C(r, 0) and the time 
dependent correlation function C+(r) = (7(0, r), which 
is essentially a Green's function. 

The spatial correlation function C x (r) can be expressed 
in terms of the link variables asQ 



C x (r) = 




(20) 
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where "path" is a straight line connecting two points a 
distance r apart in the x-direction with r fixed. Since the 
classical action (|J) is invariant under the transformation 
J x ~J X it follows that C x {r) = C x (—r). 

For the time dependent correlation function a similar 
expression can be derived. However, due to the lack of 
general particle- hole symmetry of Hamiltonian (Q), the 
correlation functions for particles C+ (r) and holes C_ (r) 
will generally be different, except for the case when /i = 
n/2 with n integer, since then statistical particle-hole 
symmetry is restored. For the single site time dependent 
correlation function one obtains 

V'epath ' 

(21) 

where " path" is a straight line connecting two points with 
fixed r a distance r apart in the (imaginary) time direc- 
tion. The site and disorder averaged correlation function 
C+(r) is then defined as 



C+(r) = [CUr)] s 



(22) 



A similar expression can be derived for C_ (r) . The single 
site function Cl(r) corresponds to the local imaginary 
time Green's function (<J>i(r)$+(0)) in the original Bose 
Hubbard model (0) . n 

According to Ref.B a scaling form for C(r, r) can be 
derived under the assumption that the long-distance and 
long-time behaviour depends on the ratios r/£ and r/£ r . 
This yields 



C(r,r)=r-^- 2 +^f(r/^r/n, 



(23) 



with a scaling function / and critical exponent r/. It 
follows that at criticality 



C x (r) 



and C+(t) ~ T- yr 



where y x and y T are given by 

y x = d—2 + z + r] 

and 

y T = [d - 2 + z + rj)/z . 



(24) 
(25) 
(26) 



Thus if one is able to determine y x and y T in the simula- 
tions, z and rj can be calculated from 



and 



z = Vx/Vr 



V = Vx - d + 2 - y x /y T . 



(27) 



(28) 



IV. MONTE CARLO ALGORITHM 

The Monte Carlo algorithm has to account for the 
zero-divergence condition of the link variables. Following 
Rcf.El we choose our update procedure to consist of lo- 
cal and global moves. In a local move four link variables 
on an elementary plaquette are changed simultaneously. 
Two of the link variables are incremented and the two 
other are decremented by one, thus fulfilling the zero di- 
vergence condition. The inverse of such a local move is 
also used, i.e. a configuration where the pluses and mi- 
nuses are exchanged. Global moves consist of changing a 
whole line of link variables extending all through the sys- 
tem by ±1 and are performed in all three directions x,y 
and r. We apply periodic boundary conditions, so the 
zero divergence condition is also fulfilled. A global move 
in the r direction amounts to injecting or destructing a 
boson, since the r component JJ of the link variables 
measures the particle number. Thus when one explic- 
itly wis hes t o keep the number of particles constant as 
in Sec. VB, global moves in the r direction have to be 



either excluded or they have to be performed pairwise, 
i.e. one creates and destructs a boson at two different 
sites simultaneously. 

For a given configuration of the link variables we then 
perform a local or global move and calculate the energy 
difference AE between these two states. We use the heat- 
bath algorithm, where the new configuration is accepted 
with probability 



w(AE) 



1 



l + exp(AE/K) ' 



(29) 



Time is measured in Monte Carlo steps (MCS), where 
one step corresponds to a sweep of local and global moves 
through the whole lattice. 

Since we are interested in equilibrium properties of 
the system equilibration is an important issue. In or- 
der to have a criterion for equilibration we ppaceed as in 
Ref.LJ where a method used for spin glassesc3 has been 
adapted. We consider two replicas a and of a system 
with the same realization of disorder and calculate the 
"Hamming" -distance, which is defined as 



t=t (i,T) 



(30) 



where to is a suitably chosen equilibration time and v — 
x,y,r. We also calculate the "Hamming" -distance for 
one replica given by 



zt 



*o (i,t) 



(31) 



When to is chosen sufficiently large, the system is equi- 
librated and one should have h^p(to) — h^(t ). Thus 
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by choosing a sequence of increasing values of to — 
10, 30, 100, 300, . . . one can give an estimate for the equi- 
libration time from the coincidence of two values. 

The equilibration time increases drastically with sys- 
tem size. We found that for smaller systems (e.g. 6x6x9) 
w 5000 MCS were sufficient, whereas the larger lattice 
sizes (10x10x25) took about 10 5 MCS. The equilibration 
time mainly depends on the system size L in the spa- 
tial directions, since the disorder, which only enters the 
transition probability of global moves in the r direction, 
enhances the acceptance rate for this type of moves. This 
made it possible to study systems with large L T such as 
8x8x200, which became necessary to avoid finite-size ef- 
fects when calculating the probability distribution of the 
local susceptibility in Sec. v[ 

One also has to consider the disorder average [. . .] av . 
The number of samples necessary to get decent statis- 
tics strongly depends on the quantity of interest. For the 
stiffness (see Section VI) some 100 samples were sufficient 
while the compressibility has much larger sample to sam- 
ple fluctuations and made it necessary to average over up 
to 3000 samples. This number of samples was also needed 
for the probability distribution of the susceptibility. 

The simulations were performed on two massively par- 
allel computers, a Parsytech GCel Transputer Cluster 
with 1024 T-805 transputer nodes and a Paragon X/PS 
10 with 140 i860 processors. 



V. BOSE GLASS PHASE 

In this section we systematically investigate the Bose 
glass phase and work out the similarities to the quan- 
tum Griffiths phase observed in random transverse Ising 
systems. We choose the "strongest" possible disorder 
A = 0.5, so that the Mott lobes totally disappear and 
only the Bose glass- and the superfluid phase remain. 
Furthermore throughout this section we assume (i = 1/2, 
which corresponds to half filling. 



A. Green's functions and the local susceptibility 

In Ref.i it is argued that due to a non- vanishing density 
of states at zero energy the disorder averaged imaginary 
time Green's function C+(t) decays algebraically 



C+(r) 



(32) 



throughout the Bose glass phase. As a consequence, the 
superfluid susceptibility x = J C(r)rfr is supposed to di- 
verge in the BG phase and not only a criticality. This 
behaviour is similar to the quantum Griffiths phase ob- 
served in other models, such as the random transverse 
Ising chaipllpr the Ising spin glass in a transverse mag- 
netic fieldliiaO, where various susceptibilities also diverge 
within part of the disordered phase. 



To verify the above scenario we calculate the imaginary 
time Green's function and the probability distribution of 
the local susceptibility in the simulations. 
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FIG. 3. The imaginary time Green's function C+(r) for 
K = 0.17 (top) and K = 0.22 (bottom) in the Bose glass 
phase for different system sizes. The dashed lines are fits 
to C+(t) = a(r VT + (L T — t) Vt ). The exponent y T is inde- 
pendent of the system size and independent of K. We get 
y T = -1.08 ±0.03. 

For the chosen disorder strength A = 0.5, the phase 
diagram of Hamiltonian (^) consists of a Bose glass phase 
for small hopping amplitudes t (i.e. small values of the 
coupling K in the classical model (ph) and a superfluid 
phase for large values of t (see Fig. g) . In order to study 
the Bose glass phase, we have to locate the critical point 
for the BG-SF transition (Point #1 in Fig. |) for the 
chemical potential (i — 0.5 we are considering. This can 
be accomplished by a finite-si ze sca ling analysis of the 
stiffness as pointed out in Sec. [II A 
coupling Kc(A = 0.5, fi 



0.5) 



We find the critical 
= 0.247 ± 0.003. Thk 
phase transition was also intensively investigated in Rcf.LJ 
in the same representation, and the value of Kc reported 
there perfectly agrees with the one we found. 

We now investigate the BG phase for couplings K < 
Kc = 0.247. We start by calculating the Green's func- 







tion C+(r) given by eq. (p2|), which is shown in Fig. g 
for couplings if = 0.17 and K = 0.22. 

One observes that C + (r) indeed decays algebraically 

as 



C+(r) 



(33) 



within the error bars and that there is no dependence 
of the exponent on the system size. Similar plots can be 
made for other values of the coupling K in the BG phase, 
and one obtains the same exponent. 

Relation ( |33| ) leads to a diverging superfluid suscepti- 
bility x ~ /C( T )^ T - A convenient quantity to further 
characterize the susceptibilty is the probability distribu- 
tion of the local, i.e. single site susceptibility. For a di- 
verging superfluid susceptibility a broad distribution of 
local susceptibilities is expected. We calculate the local 
susceptibility Xi an d the corresponding probability dis- 
tribution P(xi), where Xi is defined as 



Xi = ^c;(r) 

T = l 



(34) 



Since we expect the distribution P(xi) to be very broad, 
it is convenient to work on a logarithmic scale, i.e. we 
consider P(hiXi) rather than P(xi)- 
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FIG. 4. The probability distribution P(ln xioc) of the local 
susceptibility for various values of K < Kc = 0.247 in the 
Bose glass phase (A = 0.5, \i = 0.5). The system size is L — 6 
and L T = 200. For K = 0.19, also data for L = 4 and L = 10 
are shown, which is indistinguishable from L = 6. The dashed 
line has slope —1. 

Fig. [| shows P(lnXi) for different couplings X and 
system sizes in the Bose glass phase. One clearly ob- 
serves asymptotically (x » 1) a linear relation between 
lnP(lnx) and lnx corresponding to an algebraic decay 
P(lnx) ~ X with £ = 1. Since P(lnx) = X-P(x) this 
is equivalent to 

P( X ) ~ X- 2 • (35) 

Following the argument put forward in Ref.0 for the 
corresponding quantity in the quantum Griffiths phase of 



the transverse Ising spin glass one can thus introduce a 
dynamical exponent z: The large values of Xi originate in 
local clusters containing the site i that would have a small 
energy gap if isolated from the rest of the sample. In the 
extreme case of vanishing coupling K they simply come 
from the sites with a local chemical potential that lends a 
small energy difference between single or double occupied 
sites rii = 1 or 2, i.e. i>; = 0. Thus the probability for 
such a value of Xi is proportional to the volume of the 
system P(x) ~ L • On the other hand, one can relate 
the length scale L with the inverse of the energy scale x 
in the usual way introducing a dynamical exponent z via 



P(X)~P(L'/ X ) 



With 



we obtain 



xP(x) 



L d x -< = (Ly x ) d/Z 



lnP(lnx) 



- — lnx + const. 

z 



(36) 



(37) 



(38) 



with z = d/Q. With £ = 1 it follows that z = d = 2 
throughout the BG phase. Note that here z = d in the 
Bose glass phase and at the critical point, although the 
two exponents have their origin in different physics: the 
dynamical exponent z used in eq. (EjsJ) has its origin in 
purely local effects, whereas the critical exponent z de- 
scribes global, collective excitations at the critical point. 



0.01 




0.001 



100 

T 

FIG. 5. C-f(r) obtained by Laplace transformation (eq. 
([39J)) of P(x) for different couplings K in the BG phase. 
The dashed line is ~ r -1 . The exponent —1 is in agreement 
with the one calculated directly from the correlation function 
C+(t) (see Fig. §. 

If we suppose an exponential decay C\ ~ exp(— r/xi), 
where Xi is an effective inverse correlation length in the r 
direction, we can calculate the averaged correlation func- 
tion C+(t) = [C^(r)] av from 



C+(r) 



P(x)exp(-r/x)dx 



(39) 
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Thus we have a nice consistency check for our data by re- 
calculating C+(r) from P(x)- Fig. || shows the resulting 
curves for C+(r) for different values of K < Kq — 0.247. 
By comparison with the dotted line, one again establishes 
that C+(r) ~ r" 1 . 

The asymptotic decay of P(x) ~ X 2 that can be 
clearly seen in Fig. || indicates that the underlying 
physics is purely local. In order to see this we consider 
the imaginary time Green's function for T = in the 
limit of zero hopping t = 0. Starting from the definition 

Gt(r) = (n?|* < (r)*+(0)|n?) (40) 
and substituting <5(t) = e TH $(0)e~ TH , one obtains 



C\{r) 



e(n?)-e(n?+l)] 



(41) 



where is the ground state occupation number at site 
i. From Hamiltonian (Q) the energy difference for the 
states with rii and rij + 1 bosons at site i for t = can 
be readily calculated. We obtain e(n.j) — e(n, + 1) = 
— (rij + 1/2 — /x — Uj). Since we are working at half filling 
fi = 1/2, and Vi E [—1/2, 1/2], it follows that the on-site 
energy is minimized by n% = for Vi < and rij = 1 for 
Vi > 0. For = the states with n$ = and rij = 1 are 
degenerate. As a result we obtain 



CUr) 



e»' T Vi<0 
e Oi-i)T > o 



(42) 



From the definition Xi (uj = 0) = J °° drC l + (r) it then 
follows 



X*(^ = 0) = 



R Wl<0 
1^7 «*>0 



(43) 



and since P(v) is uniform in the interval [—1/2, 1/2], one 
finally obtains 



(44) 



From this we see that the observed decay (|35|) of P(x) 
can be explained under the assumption that the excita- 
tions in the BG phase are fully localized and the sites 
are essentially decoupled. In particular it turns out that 
increasing the coupling K between sites does not change 
the strength of the singularity in dolf ) significantly, imply- 
ing that the dynamical exponent z does not change from 
K = (z=2) to Kc- This is in contrast to what happens 
in the otherwise equivalent quantum Griffiths phase in 
random transverse Ising modelsL3'li3, where an increase 
of the coupling between the spins leads to a continous 
change of the strength of the Griffiths singularities, i.e. 
to a continously varying dynamical exponent. 

The reason for this different behaviour of the Bose 
glass phase in the disordered boson Hubbard model and 
the Griffiths phase in the random transverse Ising model 
is the fact that the former is related to an XY-model 



whereas the latter is an Ising model. The effect of 
strongly coupled clusters on the dynamics (in imaginary 
time) is much weaker for vector spins than for Ising spins 
in analogy to the relaxational dynamics in classical spin 
models within the Griffiths phaseEM 

To clarify this point and to obtain information about 
the localization of the excitations in the present model, 
we consider the participation ratio in the following sec- 
tion. 



B. Participation ratio 

In the context of localization in disordered quantum 
mechanical systems one usually studies a quantity called 
the participation ratio, which is defined with respect to 
one-particle excitations. Given the latter in the space 
representation (e.g. Vi) one obtains via the local prob- 
abilities pi — |?/>i| 2 the information about the degree of 
localization of this one-particle excitation by consider- 
ing the ratio of moments. Since the MC algorithm we 
use provides us with information on the ground state 
exclusively (not on excitations), we devise an alterna- 
tive method: we add an extra particle to the system and 
study how it distributes its local probabilities pi among 
the lattice sites i. To this end we have to keep track of 
the local particle densities of the system with and with- 
out this extra boson - the difference is then interpreted as 
the probability density pt of the extra boson. The degree 
of localization can then be read from the participation 
ratio 



PL 



N 



(45) 



For a fully localized state pi — Sij it is pl — 0(1), 
whereas for a completely delocalized state pi — 1/N one 
has p L = O(N). 
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FIG. 6. The participation ratio pl(K) defined in eq. 
as a function of K for various system sizes (A = 0.5, p, = 0.5). 
Within the Bose glass phase (K < Kc = 0.247) pl ap- 
proaches a constant. 
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To calculate Pl in the simulations we use a replica 
method. We consider two replicas a and (3 with the same 
realization of disorder and with fixed number of particles. 
Replica a is initialized with N/2 + 1 and replica (3 with 
N/2 bosons, corresponding to the case of half filling \i — 
1/2 considered in the preceding calculations. We then 
calculate the probability pi of finding the excess boson 
at site i by 



1 



E 



(46) 



The restriction of working with a fixed number of par- 
ticles in each replica makes it necessary to specifically 
consider the global moves in the r direction in the Monte 
Carlo algorithm, since they amount to either injecting or 
destructing a boson. In a first attempt one can simply 
eliminate this type of move, but it turned out that it was 
impossible to equilibrate the systems within a reasonable 
amount of computer time this way. Alternatively, one can 
also keep the number of particles constant by injecting 
and destructing a particle at two different sites simulta- 
neously. We found that by randomly choosing two sites 
where at one site a particle is injected and at the other a 
particle is destructed, the equilibration time is drastically 
reduced, enabling us to obtain equilibrium results. 
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FIG. 7. Scaling Plot of the participation ratio pl (K) with 
v = 0.9 and y = 1.0 ± 0.1. 

Fig. ^| shows pl for different system sizes with fixed 
aspect ratio 1/4 and for different couplings K. One ob- 
serves that the participation ratio increases with K un- 
til it eventually saturates. At saturation it is indeed 
as expected p = O(N), which tells us that in the su- 
pcrfluid phase the extra boson is completely delocalized. 
For lower couplings however, the participation ratio ap- 
proaches 1 with decreasing K independent of the system 
size, which indicates that the extra boson is localized on 
a single site. We remark again that the phase transi- 
tion from the Bose glass to the superfluid phase occurs 
at K c = 0.247. 



By inspection of Fig. |6| and-motivated by a recent re- 
sult within mean-field theoryc.il, one might expect that 
the participation ratio satisfies a scaling law. As an 
ansatz for a scaling relation we make 



Pl(K) 
L 2 



L- y q(LV u 8,L T /L z ) , 



(47) 



where the scaling function q and an additional exponent 
y have been introduced. 6 = (K — Kc)/Kc is the dis- 
tance from the critical point and v the correlation length 
exponent. The values Kq = 0.247 and v = 0.9 can be 
obtained fr om a scaling analysis of the stiffness accord- 
ing to Sec. |III A and the values found are in agreement 
with the values reported in Ref.El. We emphasize that 
we work at fixed aspect ratio L T /L Z = 1/4, so that the 
scaling function q just depends on one argument and the 
only variable to be determined is y. Fig. [7] shows a scal- 
ing plot of the participation ratio. The data scales nicely 
as long as the participation ratio is not saturated. We 
obtain best scaling for y — 1.0 ± 0.1. 




FIG. 8. The participation ratio p~l(K) defined in eq. 
as a function of K for various system sizes. Within the Bose 
glass phase (K < Kc = 0.247) pl approaches a constant. 

We also considered different moments of the probabil- 
ity pi. For instance, we calculated 



PL 



(48) 



which is shown in Fig. One observes that the curves 
look very similar to the ones in Fig. ^ and one obtains 
the same information about the localization of the excita- 
tions: for K < Kc — 0.247 pl approaches 1 independent 
of the system size indicating again the localization of the 
excitations in the BG phase, whereas for larger values 
of K it is Pl ~ O(N) meaning that the excitations are 
delocalized. We found that it is easier to calculate pl 
in the simulations than pl(K) 7 as it turned out that the 
equilibration time for pl is much shorter than for pl- 
This is the reason why we have more data points for pl, 
making the curves look smoother. A scaling plot of Pl 
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according to eq. ( ft7| ) is shown in Fig. g. We find the 
same exponent y = 1.0 ± 0.1 as for the other definition 
of the participation ratio in eq. (flq) . 
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FIG. 9. Scaling Plot of the participation ratio £>l (K) with 
v = 0.9 and y = 1.0 ± 0.1. 



VI. PHASE DIAGRAM AND PHASE 
TRANSITIONS FOR WEAK DISORDER 

In this section we consider the model (|J) for disorder 
strength A = 0.2. In particular we examine the phase 
diagram at the tip of the lobe, i.e. at pb = 1 (Point #1 
in Fig. |l|), where it is not clear whether there is a direct 
Mott insulator to superconductor phase transition or if 
an intervening Bose glass phase exists. Without disorder 
this phase transition is in the universality class of the 3D 
XY-model. In their scaling theory Fisher et al. argue 
that for every disorder strength A / the Mott lobes 
are surrounded by a Bose glass phase even at the tip of 
the lobe and consequently, the insulator-superconductor 
transition always occurs from the Bose glass phase. To 
check this prediction and to determine the critical ex- 
ponents we perform a finite-size scaling analysis of the 
stiffness and the compressibility and also investigate the 
Mott insulator to Bose glass transition. Furthermore we 
calculate the probability distribution of the local suscep- 
tibility P(lnx) in the Mott insulating and the superfluid 
phase. 



A. Finite-size scaling analysis at the tip of the 
lohe(p — 1) 

The relevant quantities of interest are the stiffness p 
and the compressibility k, since they allow for a distinc- 
tion of the different phases. We remark again that the 
stiffness vanishes in the Mott insulating and the Bose 
glass phase but is finite in the superfluid phase, whereas 



the compressibility is zero in the Mott insulating phase 
and finite in the Bose glass and the superfluid phase. 
Since for a finite-size scaling analysis according to Sec. 



Ill the aspect ratio L T /L Z has to be kept constant, and 
the dynamical exponent z is not known a priori, we first 
have to make a choice for z. However, by calculating 
the correlation functions ( p0| ) and ( |22] ) at criticality, we 
obtain z independent of the aspect ratio, and thus we 
have a consistency check for the assumed value of z. 
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FIG. 10. The stiffness p(0)L at the tip of the lobe 
(A = 0.2, p — 1.0). The aspect ratio is constant 
for z = 1. The intersection point is seen to be at 
K C (A = 0.2, p = 1) = 0.325 ± 0.003. 

We proceed by first assuming z = 1 which is the value 
for the pure case and remark that the following calcula- 
tion of the correlation functions shows that this is indeed 
the correct value. We calculate the stiffness p as given 
by eq. (|l l|) and perform a finite-size scaling plot guided 



by eq. 



i.e. we plot L 



d+z-2 



p vs. K for different 



system sizes but with constant aspect ratio L T /L Z = 1. 
At the critical point Lp should be independent of L, 
so at criticality we expect an intersection of the differ- 
ent curves. Fig. [l(] shows the corresponding plot. One 
clearly observes an intersection point at the critical cou- 
pling K C {A = 0.2, n = 1) = 0.325 ± 0.002. For K < K c 
the stiffness scales to zero with increasing system si 
which indicates an insulating phase. According to Ref. 
the Bose glass to superfluid transition always has z — 2, 
which suggests that the transition from the insulating 
side is from the Mott insulator rather than the Bose glass 
phase. The critical point is marked as point # 1 in Fig. 

The critical exponent v can be obtained from a scaling 
plot of the stiffness. In Fig. [ll] we plot L d+Z ~ 2 p vs. 
L x l v 8 and adjust v as to achieve best data collapse where 
S = {K - K c )jK c with K c = 0.325. One observes that 
the data scale nicely and we obtain v = 0.7 ± 0.1. 
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FIG. 11. Scaling plot of the stiffness at the tip of the lobe 
(A = 0.2, fi = 1) with K c = 0.325 and v = 0.7 ± 0.1. 

We proceed by calculating the compressibility k ac- 
cording to eq. (|l8|). Fig. 12 shows a finite-size scaling 
plot of k with z - 
point close to K 



1. One again observes an intersection 
0.325, even though the intersection 
point is not as good as the one for the stiffness. This 
is due to fact that the compressibility has much larger 
sample to sample fluctuations than the stiffness and so 
it is more difficult to obtain good statistics. Note that 
the number of samples for the compressibility already is 
2500, whereas for the stiffness 250 samples were suffi- 
cient. The large fluctuations also limited the range of 
system sizes we were able to simulate for the compress- 
ibility, so we only have data up to systems of size 8x8x8. 
For couplings K < Kc — 0.325 we see that Lk decreases 
with increasing system size, which indicates that in the 
thermodynamic limit K = Q. 
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FIG. 12. The compressibilty at the tip of the lobe 
(A = 0.2, n — 1.0) for different sytem sizes. The aspect ratio 
is constant for z = 1. The intersection point coincides with 
the one for the stiffness in Fig. [l^ within the error bars. 

Considering the above results for the stiffness and the 
compressibility, we conclude that below Kc — 0.325 we 
have indeed a Mott insulating phase and for K > Kc a 



superconducting phase, in particular there is no sign of 
an intervening Bose glass phase. 
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FIG. 13. The correlation functions C x (r) (top) and C+(r) 
(bottom) for system sizes 8x8x16 and 10x10x10 at the crit- 
ical point Kg = 0.325. The dashed lines are fits to 
C x (r) = a(r y " + (L — r) y °>) and C+(r) = b(r VT + {L T - r)^) 
respectively. We obtain z = yx/t/r ~ 1 independent of the 
shape of the sample. 

We can now go for the consistency check of our as- 
sumed value of z = 1. Therefore we measure the spatial- 
and imaginary time correlation functions as given by eq. 
( pp| ) and (p2|). From these functions we obtain the expo- 
nents y x and y T and then z is simply given by eq. (|27j), 
i.e. z — Hx/Ut- Fig. [Dj] shows the correlation functions 
for different system sizes at criticality Kc — 0.325. By 
fitting to the functions C + (t) — o,(t~ Vt + (L T — t)~ Vt ) 
and C x (r) = a'(r~ V:c + (L — r)~ V:c ) respectively we ob- 
tain y T = 1.11 ± 0.02, y x = 1.07 ± 0.02 for system size 
10x10x10 and y T = 1.10 ± 0.02, y x = 1.10 ± 0.02 for 
system size 8x8x16. This yields z — y x /y T = 1 in both 
cases, confirming the value used in the scaling analysis of 
the stiffness and the compressibility and demonstrating 
that z is indeed independent of the shape of the sam- 
ple. Additionally from y x and y T the critical exponent 
r] can be obtained according to eq. (|2q), which yields 
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7/ = 0.1 ±0.1. 

To summarize the above scaling analysis, we find a 
direct Mott insulator to superconductor transition with 
no intervening Bose glass phase at the multicritical point 
fi = 1. The critical exponents are z = 1, v = 0.7 ± 0.1 
and r\ = 0.1 ± 0.1. These are very close to the exponents 
obtained fox,the 3D XY-model both analytically from e- 
expansionsB v = 0.669 ±JX002, r\ = 0.033 ±0.004 and by 
Monte Carlo simulationsu where v = 0.667 and r\ w 0.02 
were obtained. Thus we conclude that the universality 
class of the transition is not changed by weak disorder. 
However, this might be different for stronger disorder, 
and for instance for A = 0.4 we get an effective exponent 
z = 1.4. 

Furthermore we want to comment on the result v = 
0.7. ThiSj-sialue of v seems to contradict the inequality 
v > 2/c&3'c3, where d is the dimension in which the sys- 
tem is disordered (i.e. here d = 2 ), which should be 
applicable to a broad class of phase transitions that can 
be triggered by the strength of the disorder. However, by 
comparing the critical couplings for the MI-SF transition 
for the pure case K C (A = 0, fi = 1) = 0.333 ± 0.003 and 
weak disorder K C {A = 0.2, (i = 1) = 0.325 ± 0.002 it 
becomes evident, that the position of the tip of the lobe 
only depends very weakly on the disorder strength or is 
even independent of it, at least for small disorder. Thus 
this transition cannot be triggered by varying the disor- 
der strength and therefore the inequality possibly does 
not apply here. We also want to mention the possibil- 
ity, that our data are not yet in the asymptotic scaling 
regime and therefore the measured exponent v is only an 
effective exponent for small length scales, although we 
find this very unlikely due to the quality of the scaling 
plots. 

We remark that we also find an intersection point in 
the finite size scaling plot of the stiffness at K = 0.32 
(which is close to the critical coupling Kc = 0.325 we 
find for z = l) when assuming z = 2, even though this in- 
tersection point is not as convincing as the one for z = 1. 
However, as pointed out above, by calculating the corre- 
lation functions we always obtain z = 1, independent of 
the shape of the sample. Furthermore there is no inter- 
section point in the plot of the compressibility for z = 2. 
Thus a transition with z = 2, which could indicate a 
transition from the Bose glass phase, is ruled out by the 
above results. 



B. The Mott Insulator - Bose glass transition 

As we have seen the BG phase is quite similar to the 
quantum Griffith phase occuring in random transverse 
Ising models. One can extend the analogy to the other 
phases, too, and formulate the following dictionary (see 
Table |): 

One can identify the Mott insulating and paramagnetic 
phases as the disordered phases, the Bose glass and the 



Griffiths phase as the weakly or locally ordered phases 
and the superfluid and ferromagnetic phase as the or- 
dered phases. The transition from the insulating BG 
phase to the superconducting phase as well as the transi- 
tion from the paramagnetic Griffiths phase to the ferro- 
magnetic phase are both second order, whereas the MI- 
BG and PM-Griffiths "transitions" are in fact not real 
phase transitions: only some susceptibility starts to di- 
verge due to local effects, which for the Bose glass phase 
also implies that it becomes compressible. As in the con- 
text of the Griffiths phase, where a continuously increas- 
ing (coming from the disordered phase) dynamical expo- 
nent leads to this divergence, no diverging length scale is 
connected with the closing of the gap upon entering the 
Bose glass phase (coming from the Mott insulator). 
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FIG. 14. Number of particles per site vs. chemical po- 
tential at K = 0.19 for different system sizes and different 
disorder strengths A. In the Mott phase the particle number 
is fixed at an integer value. For increasing disorder (A = 0.4) 
the Mott phase shrinks. 

As a consequence we cannot locate the MI-BG transi- 
tion via finite size scaling. To get an idea of the phase 
diagram, we measure the particle number rij per site 
as a function of the chemical potential fi at fixed cou- 
pling K — 0.19. This corresponds to moving along 
the line through point #2 in Fig. |l|. We choose the 
value of K = 0.19 as to be below the critical coupling 
K C (A = 0.2, fi = 0.5) = 0.20 for the BG-SF transition 
(Point #3 in Fig. 0). Such a plot is shown in Fig. [l4| 
for K = 0.19. The different curves correspond to dif- 
ferent system sizes. One observes that starting at the 
chemical potential fj, = 0.65 the particle number per site 
increases with roughly constant slope. This means, that 
the compressibility k = dn/dfi is indeed finite in a cer- 
tain range of the chemical potential. For /i ks 0.73 the 
particle number per site reaches one, meaning that we 
are in the Mott insulating phase. Increasing the chemi- 
cal potential further does not change the number of par- 
ticles, which means that the compressibility vanishes as 
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expected. For higher values of /x, we leave the Mott phase 
and again have a phase with non vanishing compressibil- 
ity. For stronger disorder (A = 0.4) one observes a sim- 
ilar behaviour, but the range of the chemical potential 
where the particle number is constant is smaller. This 
indicates that with increasing disorder the Mott lobes 
indeed shrink, as the disorder reduces the energy gap in 
the MI phase. 

In Fig. |l5] we show the number of particles vs. chemi- 
cal potential near the tip of the lobe. It can be seen that 
at the critical point Kc(A = 0.2, fi = 1) = 0.325 there 
is no longer a range of the chemical potential where the 
number of particles is constant demonstrating that the 
MI phase disappears at I<c- 
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FIG. 15. Number of particles per site vs. chemical poten- 
tial for A = 0.2 near the tip of the lobe for system size 6x6x18, 
demonstrating that the Mott lobe disappears at Kc = 0.325. 




FIG. 16. The compressibility k as a function of the chemi- 
cal potential fi for the Mott insulator to Bose glass transition 
at fixed K = 0.19 (A = 0.2). 



at the MI-BG transition (which does not diverge here, as 
we mentioned above) is smaller than the system size L, 
the compressibility k will not depend significantly on L. 
Indeed we observe a much stronger dependence on L t au, 
which is due to the closing of the gap at the MI-BG tran- 
sition. In Fig. 16 we demonstrate this effect by showing 



data for increasing system size in the imaginary time di- 
rection and constant system size in the space direction. 
The curves intersect at yi = 0.73, which corresponds to 
the left edge of the Mi-plateau in Fig. [lj, i.e. the MI-BG 
transition point. For increasing L T the curves become 
steeper and one might guess that in the limit L T — » 00 
the compressibility jumps discontinously from k = for 
/i > 0.73 (in the Mi-phase) to a finite value (k w 0.5) for 
/i < 0.73 (in the BG phase). Unfortunately the sample 
to sample fluctuations of k are extremely large so that 
we could not get good enough statistics for larger values 
of L T . 

A discontinuity in the compressibility would support 
the conjecture that the MI-BG transition is first-ordem. 
At a generic first order phase transition one expects a 
phase coexistence, which in our case means the coex- 
istence of the compressible BG and the incompressible 
MI phase at fi «r-Q,.73. Such a possibility can be stud- 
ied systematicallyEa via the calculation of the probability 
distribution P(k) of the compressibility in finite systems. 
If the MI-BG transition is first order and shows phase co- 
existence, one expects a double peak structure in P(k) 
with a peak at zero (from MI phase) and at ss 0.5 (from 
the BG phase). However, our data for P(n) (see Fig. |l7j ) 
show only a single peak moving from n — to k ~ 0.5 
when increasing /i from below to above the transition, 
which excludes phase coexistence. 
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FIG. 17. The probability distribution P(k) of the com- 
pressibility for system size 6x6x9 at /1 = 0.73 and K — 0.19 
(A = 0.2). Note the logarithmic scale of the y-axis and the 
exponential tail of the distribution. 



We also calculate the compressibility directly as func- 
tion of )jl at the same value of K = 0.19 as above, which 
is shown in Fig. [lfj. As long as the typical length scale £ 
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C. The local susceptibility in the Mott insulating 
and superfluid phase 

Finally we consider the probability distribution P(ln x) 
of the local susceptibility in the Mott insulating and the 
superfluid phase. We fix the chemical potential at (J, = 1, 
i.e. we are at the tip of the lobe and calculate P(ln%) for 
different couplings K. The resulting curves are shown 
in Fig. From the preceding section we know that 

the Mott insulator to superfluid transition for A = 0.2 
and fj, = 1.0 is at Kc — 0.325. For the smallest value 
of K = 0.250, which clearly is in the Mott insulating 
phase, one observes that P(lnx) has no tail as in the 
Bose glass phase but rather abruptly decays to a very 
small value. This indicates that in the MI phase one 
has indeed a finite energy gap. At criticality, i.e. at 
K = 0.325, the imaginary time Green's function decays 
algebraically with an exponent close to one (« 1.1, see 
Fig. |l3|). Since J drC(r) ~ J dxxP(x)i we again expect 
a probability distribution P(lnx) with a 1/x tail, which 
is verified by comparison with the dashed line in Fig. [l^. 
For higher values of K in the superfluid phase the curves 
are shifted to the right, and one again has an algebraic 
tail. 



K=0.250 
K=0.325 
K=0.350 
K=0.500 
K=0.700 
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x=lnxi 0c 

FIG. 18. The probability distribution P(ln Xioc) of the local 
susceptibility in the Mott insulating and the superfluid phase 
as a function of K for A = 0.2 and fj, = 1. The MI-SF 
transition is at K c = 0.325. In the MI phase (K = 0.250) 
the distribution P(lnx) has no tail, indicating a finite energy 
gap. At criticality and in the SF phase one recovers a tail. 
The dashed line has slope -0.9. 



The occuring singularities can be parametrized by a dy- 
namical exponent z — 2 which is constant in the Bose 
glass phase. This dynamical exponent z is equal to the 
critical dynamical exponent z — d — 2 for the Bose glass 
superfluid transition and it is an open question whether 
this equality holds in general, since the exponents have 
their origin in different physics. We remark that this 
equality was-also found to hold in the random transverse 
Ising-chainO and in the random transverse Ising spin 
glassEl 

By calculating the participation ratio we were able to 
show that the excitations in the Bose glass phase are fully 
localized. This is in contrast to the random transverse 
Ising systems, since in these models rare strongly coupled 
clusters of spins act collectively giving rise to a divergent 
susceptibility. Furthermore we find the participation to 
obey a scaling relation. 

The investigation of the phase transition at the tip of 
the lobe for integer boson densities showed a direct Mott 
insulator to superconductor transition for weak disorder. 
The critical exponents we obtain are z = 1, v = 0.7 ±0.1 
and rj = 0.1 ± 0.1 which agree with the exponents of the 
pure 3D XY model, so the universality class of the tran- 
sition is not altered by weak disorder. This conclusion 
is supported by recent mean field renormalization group 
studies and scaling arguments put forward mE3. 
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VII. SUMMARY AND DISCUSSION 

Concluding we have shown that the Bose glass phase 
shares many features with the Griffiths phase in random 
transverse field Ising systems. We find an algebraic decay 
of the imaginary time Green's function, a broad proba- 
bility distribution of the local susceptibility with an 1/x 2 
tail and thus a diverging global superfluid susceptibility. 



TABLE I. Corresponding phases in the disordered boson 
Hubbard model (d. B-H-M) and the random transverse Ising 
chain (RTIC). p is the superfluid stiffness, m the spontaneous 
magnetization and xioc the local susceptibility. 
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